Relatorio de obtenção dos resultados

library(gtools)
library(doMC)
library(GUILDS)
library(lme4)
library(merTools)
library(magrittr)
library(gridExtra)
library(ggplot2)
library(stringr)
library(tidyr)
library(plyr)
library(purrr)
library(dplyr)

Seleção e Preparo dos Sítios de Amostragem

Dados Disponíveis

Seleção dos trabalhos fitossociológicos

Filtros gerais:
i) effort >= 1ha; ii) DBH>=5cm; iii) ano dos dados >= 2000;

Filtros condicionais:
i) state %in% Rio de Janeira, Rio Grande do Sul -> ano dos dados >= 1990
ii) state %in% Bahia, Goiás, Mato Grosso do Sul -> ano dos dados >= 2000
iii) para as demais regiões -> ano dos dados>= 1995
iv) Exceções a esse esquema ocorreram quando trabalhos foram feitos antes do ano de 2000 em grandes áreas de regiões protegidas (>1000ha) ou em antigos campi universitários.

### dados
df_references <- read.csv(file="/home/danilo/Documentos/Doutorado/artigo_mestrado/1A.P_MOVER/dados_brutos/references - TreeCo.csv",
                          as.is = TRUE,header =T, na.strings = c("","NA"))
### filtro geral
df_1ofiltro <- df_references %>% filter(method=="parcelas" &
                                         grepl("*contiguous*",arrangement) & 
                                        effort_ha>=1 &
                                         grepl("*yes*",status) & 
                                         grepl("ok*",status_diagnostico) &
                                         grepl("Atlantic_Forest*",domain) &
                                        dbh_cutoff %in% c("PBH>=15.0cm","PBH>=15.7cm",
                                                          "DBH>=5.0cm","DGH>=5.0cm",
                                                          "DBH>=4.8cm",
                                                          "DBH>=5.0cm&H>300cm","DBH>=5.0cm&H>500cm", 
                                                          "DGH30>=5.0cm") )
### filtro universidades e campi
df_simulacao_s_UCprotecao.integral <- df_1ofiltro %>% 
  filter(UC_area_ha >= 1000 & 
         Unidade_de_conservacao == "universities and research centers")
### filtro estados
df_simulacao_filtro.estate <- filter(df_1ofiltro,state %in% c("RJ","RS") & year >= 1990)
df_simulacao_filtro.estate %<>% rbind(.,filter(df_1ofiltro,state %in% c("BA","GO","MS") & year >= 2000)) %>% 
  unique
df_simulacao_filtro.estate %<>% rbind(.,filter(df_1ofiltro,!(state %in% c("BA","GO","MS","RJ","RS")) & year >= 1995)) %>% 
  unique
### unindo e salvando
df_ref.S_UCprotecaoIntegral <- rbind(df_simulacao_s_UCprotecao.integral,df_simulacao_filtro.estate) %>% 
  distinct()

SADs obs

• A SAD obs conta com a abundância dos indivíduos mortos em pé, que foram desconsiderados para a parametrização das simulações.

# SAD obs
df_SAD.obs <- read.csv(file="/home/danilo/Documentos/Doutorado/artigo_mestrado/1A.P_MOVER/dados_brutos/abundances.csv",header = TRUE,as.is = TRUE)
df_SAD.obs$SiteCode %<>% as.factor()
df_SAD.obs %<>% filter(species.correct != "Mortas") %>% select(SiteCode,RefID,ordem,species.correct,N) #%>% group_by(SiteCode,RefID,ordem) %>% nest
names(df_SAD.obs)[2] <- "refID"
# merge
df_dados.disponiveis %<>% inner_join(x=.,
                                     y=ddply(df_SAD.obs,"SiteCode",summarise,
                                            Ntotal=sum(N), Stotal=length(species.correct)),
                                     by="SiteCode")

Raster de Paisagem

• Os rasters de paisagem atualizados são para áreas > 5km2.
  
• A resolução dos arquivos é de 30x30m e o lado da paisagem é de 800 pixels.
  
• Selecionei os 200 pixels centrais, assim obtive paisagens de 6x6km.

func_tif.png <- function(file){
  raster <- raster(file) #leitura do raster
  # raster <- raster(df_dados.disponiveis$tif.name[df_dados.disponiveis$SiteCode=="MGuberl3"])
  mat_raster <- matrix(data = getValues(raster)/100) #convertendo para matrix
  dim(mat_raster) <- dim(raster)[1:2]
  #recorte de 6x6km: os 200x200 pixeis centrais
  mat_raster <- mat_raster[301:500,301:500]
  squash::savemat(x = mat_raster, filename = gsub(".tif",".png", file)) #salvando como png
}
registerDoMC(3)
# df_dados.disponiveis %>% filter(is.na(tif.name))
a_ply(df_dados.disponiveis$tif.name,1,func_tif.png,.parallel = TRUE) #gera 113 .png, oriundos dos .tif

• Ajustei a resolução das imagens de tal modo que o total de pixels na paisagem = DA_{obs} * Area_{paisagem}

func_png.ajust <- function(file, densidade, A_landscape=3600){ # atualizar para o pacote 'magick'
  system(paste(
    "convert ",file, " -resize ", densidade*A_landscape,"@ ", file,  
    sep = ""
  ))
}
df_dados.disponiveis %<>% mutate(DA=Ntotal/effort_ha) # portanto exclui os indivíduos mortos
registerDoMC(3)
a_ply(df_dados.disponiveis,1,function(x) func_png.ajust(file = x$png.name, densidade = x$DA), .parallel = TRUE) # ajuste da resolucao em funcao da densidade

• então aplico o seguinte tratamento nas imagens:

1. pixels com cobertura vegetal >= 70% em unidades de habitat e abaixo unidades de não-habitat
   
2. função focal para tapar buracos e apagar fragmentos com apenas 1 unidade de habitat
   
3. função ‘area simulada’, que marca unidades de habitat que serão seguidas na simulação coalescente:
   

1. define-se uma janela de observação; A_{janela de observação} ~ 2.25 A_{sitio de amostragem} ;
   
2. esta janela é posicionada na região central da paisagem ;
   
3. a função marca as unidades de habitat pressupondo uma espiral quadrada divergente que se inicia no pixel central da paisagem ;
   
4. unidades de não-habitat não são marcadas.

f_area.simulada <- function(matriz, N){
#@ matriz :: objeto matriz que representa
#@ N:: tamanho da amostra de indivíduos
  # Janela de observação
  d <- ceiling(sqrt(N)*(3/4)) # metade do lado do janela de observação
  l <- ceiling(dim(matriz)[1]/2) # linha central da paisagem
  c <- ceiling(dim(matriz)[2]/2) # coluna central da paisagem
  # define uma janela central na paisagem onde o for sera aplicado
  m_temp <- matriz[(l-d):(l+d),(c-d):(c+d)]
  if(length(m_temp[m_temp==1]) < N) { 
    stop("habitat insuficiente na janela de observação")
  # tambem deve estar errado, pois janela de observacao ~ 2.25A_sitio
  } else if (length(m_temp[m_temp==1]) == N) { 
    stop("area amostral igual janela de observacao")
  # paisagens que devem estar adequadas para os métodos:
  } else { 
    # posição de cada elemento da janela de observação, por coluna em ordem crescente
    col_cresc <- which(m_temp==m_temp, arr.ind = T)
    # idem na ordem contrária
    col_decre <- col_cresc[dim(col_cresc)[1]:1,]
    # posicao de cada elemento, por linha em ordem crescente
    row_cresc <- col_cresc[order(col_cresc[,1],decreasing = FALSE),] 
    # idem ao controario
    row_decre <- row_cresc[dim(col_cresc)[1]:1,] 
    # (nrow - 1)/2; exclui a posição dos elementos da coluna central da janela de observação
    ciclo <- (dim(m_temp)[1]-1)/2 # 
    l_mat_index <- list() #lista que vou usar dentro do for
    dim_temp <- dim(m_temp)[1]
    for(i in 1:ciclo){
      a1 <- col_cresc[col_cresc[,"col"]==i,] #considerando o primeiro ciclo: 1a coluna em ordem crescente
      b1 <- row_cresc[row_cresc[,"row"]==dim_temp+1-i,] #última linha em ordem crescente 
      c1 <- col_decre[col_decre[,"col"]==dim_temp+1-i,] #última coluna em ordem reversa
      d1 <- row_decre[row_decre[,"row"]==i,] #primeira linha em ordem reversa; os demais ciclos são com a segunda coluna, penúltima linha, penúltima coluna e segunda linha, etc  
      l_mat_index[[i]] <- do.call(rbind,list(a1,b1,c1,d1)) #ao final de cada ciclo eu concateno tudo em uma única matriz
    }
    l_mat_index[[(dim_temp+1)/2]] <- col_cresc[col_cresc[,"col"]==(dim_temp+1)/2,] #a coluna central deve ser a última
    mat_ref <- unique(do.call(rbind, l_mat_index)) #remocao de repeticao. Obtenho uma sequência de elementos que descreve uma espiral quadrada convergente
    length_ref <- length(m_temp[mat_ref][m_temp[mat_ref]==1]) #variável para indexação:
    m_temp[mat_ref][m_temp[mat_ref]==1][(1+length_ref-N):length_ref] <- 2 #os N últimos elementos que são iguais a 1 e troco por 2
    matriz[(l-d):(l+d),(c-d):(c+d)] <- m_temp #substituo a matriz de volta
    return(matriz)
  }
}

f_mat.tri <- function(png, abund){ #png.file, número de indivíduos presente 
  janela <- matrix(1,3,3) 
  raster <- raster(png)
  mat <- matrix(getValues(raster)/255, ncol = ncol(raster), nrow = nrow(raster))
  raster_binario <- raster( matrix(nrow = nrow(mat), ncol = ncol(mat), sapply(mat, function(x) ifelse(x >= 0.7, 1, 0)) ) ) 
  func_focal <- function(x) ifelse(sum(x[x==1]) >= 5, 1, x[5])
  binario.focal <- as.matrix( focal(raster_binario, janela, func_focal, pad=TRUE, padvalues = 0))
  mat_tri <- try(f_area.simulada(matriz = binario.focal, N = abund))
  if(class(mat_tri) == "matrix"){ 
    try(write.table(x = mat_tri, 
                    file = gsub(".png",".txt", png),
                    sep = " ", row.names = FALSE, col.names = FALSE))
  }
}
df_dados.disponiveis$png.name %<>% as.character
registerDoMC(3)
a_ply(df_dados.disponiveis,1,function(X) f_mat.tri(png = X$png.name,abund = X$Ntotal),.parallel = TRUE)

• Essa rotina gera paisagens onde a comunidade local é representada exatamente por um quadrado (figura 1).
  
• Em outros casos a particular configuração espacial não permite construir a comunidade local como uma região quadrada (figura 1).

Figura 1 Na esquerda consideramos a comunidade local como quadrada; na direita consideramos a comunidade local aproximadamente como um quadrado

• Consideramos que os métodos são uma boa aproximação quando todo o habitat da comunidade local esta conectado.
• Para as paisagens onde o habitat da comunidade local esta fragmentado em fragmentos florestais vizinhos, eu removi o habitat isolado e marquei quantidade equivalante de habitat adjacente à comunidade local:

par(mar=c(0.2,0.2,1.2,0.2))
i <- 1 #nrow = 9
N <- df_conserto.land$N[i]
paisagem <- as.matrix(read.table(df_conserto.land$txt.name[i],header = TRUE,sep = " ",as.is = FALSE))
# janela de observação:
d <- ceiling(sqrt(N)*(3/4))  # metade do lado do janela de observação
l <- ceiling(dim(paisagem)[1]/2) # linha central da paisagem
c <- ceiling(dim(paisagem)[2]/2) # coluna central da paisagem
# janela de observação:
janela_observao <- paisagem[(l-d):(l+d),(c-d):(c+d)]
###############################
########### região para remover
###############################
## objetos
#visualização
image(janela_observao,main=df_conserto.land$SiteCode[i],col=terrain.colors(12,rev = TRUE))
coordenada_ <- locator(1) %>% unlist()
l_ <- round(coordenada_[1]*nrow(janela_observao))
c_ <- round(coordenada_[2]*ncol(janela_observao))
n_ <- 5
## auditoria visual
image(janela_observao[(l_-n_):(l_+n_),(c_-n_):(c_+n_)],
      main=df_conserto.land$SiteCode[i],col=terrain.colors(12,rev = TRUE))
## para aqueles sem o habitat central como 1:
# janela_observao[(l_-n_):(l_+n_),(c_-n_):(c_+n_)] <- 2
## janela de conversao de 2->1
#1
p_arrumar <- janela_observao[(l_-n_):(l_+n_),(c_-n_):(c_+n_)]
rm_ <- p_arrumar[p_arrumar==2] %>% length()
p_arrumar[p_arrumar==2] <- 1
janela_observao[(l_-n_):(l_+n_),(c_-n_):(c_+n_)] <- p_arrumar
#2
# p_arrumar <- janela_observao[(l_-n_):(l_+n_+2),(c_-n_):(c_+n_)]
# rm_ <- length(p_arrumar[p_arrumar==2]) + rm_
# p_arrumar[p_arrumar==2] <- 1
# janela_observao[(l_-n_):(l_+n_+2),(c_-n_):(c_+n_)] <- p_arrumar
######################################
########### região para incluir 1 -> 2
######################################
# length(rm_)
## 
image(janela_observao)
coordenada_ <- locator(1) %>% unlist()
l_ <- round(coordenada_[1]*nrow(janela_observao))
c_ <- round(coordenada_[2]*ncol(janela_observao))
n_ <- 2
janela_observao[(l_-n_):(l_+n_),(c_-n_):(c_+n_)] # %>% image
# para aqules com ponto central = 1
# janela_observao[l_+1,c_+1] <- 1
# exemplo: o ponto de inicío é [4,3] e o ponto l_,c_ é [3,3] portanto em termos de l_,c_: 
p_inicio <- c(l_+1,c_+1)
# correr pela ___ até o ponto:
image(janela_observao)
coordenada_ <- locator(1) %>% unlist()
l_ <- round(coordenada_[1]*nrow(janela_observao))
c_ <- round(coordenada_[2]*ncol(janela_observao))
n_ <- 2
janela_observao[(l_-n_):(l_+n_),(c_-n_):(c_+n_)] # %>% image
p_final <- c(l_+1,c_)
# trocando 1 -> 2
p_rm_ <- janela_observao[p_inicio[1]:p_final[1],p_inicio[2]:p_final[2]] %>% length()
# if
p_rm_ < rm_
 janela_observao[p_inicio[1]:p_final[1],p_inicio[2]:p_final[2]] <- 2
 rm_ <- rm_ - p_rm_
# else 
 janela_observao[p_inicio[1]:p_final[1],p_inicio[2]:p_final[2]][1:rm_] <- 2

## outro
# pontos_p_marcar <- janela_observao[p_inicio[1]:p_final[1],p_inicio[2]:p_final[2]]
# pontos_p_marcar[pontos_p_marcar==1][1:length(rm_)] <- 2
# janela_observao[p_inicio[1]:p_final[1],p_inicio[2]:p_final[2]] <- pontos_p_marcar 
# image(janela_observao)
# auditoria e gravar
image(paisagem[(l-d):(l+d),(c-d):(c+d)],
      main=paste(df_conserto.land$SiteCode[i], "antes"),
      col=terrain.colors(12,rev = TRUE))
image(janela_observao,
      main=paste(df_conserto.land$SiteCode[i], "consertada"),
      col=terrain.colors(12,rev = TRUE))
paisagem[(l-d):(l+d),(c-d):(c+d)] <- janela_observao
image(paisagem,main=df_conserto.land$SiteCode[i],col=terrain.colors(12,rev = TRUE))
write.table(paisagem,
            file=df_conserto.land$txt.name[i],
            row.names = FALSE,col.names = FALSE)

Um exemplo de paisagem modificada esta na figura 2.

Figura 2 Exemplo de modificação realizada para poder aproximar a paisagem pelos métodos utilizados.

Por conta da função focal a paisagem esta cercada por NAs (linhas e colunas marginais), então removo-as utilizando:

f_NA_zero <- function(path_){
  mat_paisagem <- read.table(file=path_,header = FALSE)
  mat_paisagem[is.na(mat_paisagem)] <- 0
  write.table(x = mat_paisagem,file = path_,sep = " ", row.names = FALSE, col.names = FALSE)
}

Parametrização

A relação entre os parâmetros de diversidade pode ser obtida a partir da igualdade proposta por por Vallade & Houchmandzadeh (2003) \(\theta = \frac{U (J_M - 1)}{1-U}\). Para isso aproximamos o número de indivíduos na metacomunidade como \(J_M = A_{landscape} DA_{obs} p\).

O parâmetro m é a probabilidade de um evento de colonização na comunidade local ser por um propágulo da metacomunidade (Hubbel 2001). Uma vez que toda unidade de habitat apresenta igual probabilidade de colonização, podemos definir m como a média das probabilidades de cada sítio na comunidade local ser colonizada por um imigrante (eq 1; Chisholm & Lichstein 2009).

\[ eq 1: m = \frac{1}{A} \int \int_A m_{x,y} dxdy \]

Aproximamos as comunidades locais como áreas quadradas com lado L, que pode ser obtido por \(L = \sqrt{10000\frac{J}{DA}}\). Assim podemos reescrever a equação 1 como:

\[ eq 2.a : m = \left(\frac{1}{L} \int\limits_{-L/2}^{L/2} m_{x}(x)\mathrm{d}x \right)^2 \]

\[ eq 2.b :m_{x} = 1 - \int\limits_{-L/2}^{L/2} K(x-y) \mathrm{d}y \]

Onde K é a função de dispersão. Na simulação coalescente a dispersão resulta do sorteio indepentende de duas distribuições de Laplace em eixos ortogonais centradas no sítio vago (figura X). Se consideramos a distribuição de Laplace como \(K(x) = \frac{\alpha}{2} e^{-|\alpha x|}\) (idem para o eixo y), então:

\[eq3:m = (\frac{1-e^{-\alpha L}}{\alpha L})^2\]

Onde \(\alpha = 1/b\), b é o parâmetro escalar da distribuição de Laplace que pode ser escrito em função do desvio-padrão da distribuição de Laplace (d) \(b = d/ \sqrt{2}\). O desvio padrão corresponden à distância média de dispersão (Clark et al. 1999). Podemos reescrever a equação de m em funçao de d:

\[ eq4: m = d \frac{1 - e^{-\frac{\sqrt{2} L}{d}} }{\sqrt{2} L} \]

Com a equação 4 podemos calcular m a partir do desvio padrão da função de dispersão. Essa equação é válida para o processo de dispersão em paisagens homogêneas. Na simulação coalescente, uma vez que sorteamos um progenitor e este estaria presente em uma unidade de não habitat, o sorteio é refeito até que o progenitor esteja em uma unidade de habitat. Uma aproximação do efeito da fragmentação na simulação no m calculado a partir da equação 4 pode ser obtido por:

\[eq.5: m' = \frac{mp}{1 - (1-p)m} \]

Onde p é a porcentagem de cobertura vegetal na paisagem. Caso seja necessário calcular d a partir de m, utilizamos o ramo principal da função W de Lambert (\(W_{0}\)):

\[ eq6 : d = \frac{\sqrt{2} L m}{m W_{0}(- \frac{e^{-1/m}}{m} ) + 1} \]

Proporção de cobertura vegetal

Para calcular a proporção de cobertura vegetal utilizei a função:

f_tree.cover <- function(file_path){
  mat_paisagem <- read.table(file=file_path,sep=" ",header=TRUE)
  tree.cover <- as.vector(mat_paisagem)
  tree.cover <- tree.cover[!is.na(tree.cover)]
  p <- 1 - length(tree.cover[tree.cover==0])/length(tree.cover)
  return(p)
}
registerDoMC(4)
df_simulacao$p <- aaply(df_simulacao$txt.name,1,f_tree.cover,.parallel = TRUE)
# 

A relação entre proporção de habitat (p) e riqueza observada (S) esta na figura 3.

Figura 3 Distribuição de p, s e relação entre as duas variáveis. As linhas verticais em vermelho marcam respectivamente o quantil de 25, 50 e 75% da distribuição.

Parâmetro de dispersão

Parametrizamos por k, a proporção de propágulos que permanecem até a vizinhança imediata da planta progenitora. A estimativa da distância média de dispersão para o respectivo k para cada paisagem foi obtida pelas funções:

library(rmutil)
library(lamW)
qkernel<- function(sigma, kernel, p, density=20852/50, npoints = 1e5){
    kernel <- match.arg(kernel, choices=c("normal","gaussian","laplace","uniform"))
    d_ind_MA  <- 100/sqrt(density)
    if(kernel=="laplace"){
        b_laplace <- sigma / sqrt(2)
        X_laplace <- d_ind_MA * round(rlaplace(npoints, s=b_laplace) / d_ind_MA) 
        Y_laplace <- d_ind_MA * round(rlaplace(npoints, s=b_laplace) / d_ind_MA)
        dist_laplace <- sqrt(X_laplace^2+Y_laplace^2)
        result <- quantile(dist_laplace, p)
    }
    if(kernel=="normal"|kernel=="gaussian"){
        b_norm <- sigma 
        X_norm <- d_ind_MA * round(rnorm(npoints, sd=b_norm) / d_ind_MA)
        Y_norm <- d_ind_MA * round(rnorm(npoints, sd=b_norm) / d_ind_MA)
        dist_norm <- sqrt(X_norm^2+Y_norm^2)
        result <- quantile(dist_norm, p)
    }
    if(kernel=="uniform"){
        b_unif <- sigma/2
        X_unif <- d_ind_MA * round(runif(npoints, min = -b_unif, max = b_unif) / d_ind_MA)
        Y_unif <- d_ind_MA * round(runif(npoints, min = -b_unif, max = b_unif) / d_ind_MA)
        dist_unif <- sqrt(X_unif^2+Y_unif^2)
        result <- quantile(dist_unif, p)
    }
    return(unname(result))
}

sigkernel <- function(kernel, p, distance, density=20852/50,
                      npoints =1e5, sigma.min = 1, sigma.max= 100){
    f1 <- function(x) distance - qkernel(x, kernel, p, density, npoints)
    uniroot( f1 , lower = sigma.min, upper = sigma.max)
}

Utilizamos os percentis: c(0.99,seq(0.95,0.05,-0.05)) (figura 4). Obtenção dos percentis:

# dados
percentil <- c(0.99,seq(0.95,0.05,-0.05))
df_simulacao %<>% left_join(x=.,
                            y=expand.grid(SiteCode = df_simulacao$SiteCode, k = percentil),
                            by="SiteCode")

df_simulacao$kernel_type <- "laplace"
df_simulacao$kernel_code <- "2"
df_simulacao %<>% mutate(DA=Ntotal/effort_ha, # densidade desconsiderando as espécies mortas
                         dist_0 = 100/sqrt(DA)) # distância entre indivíduos vizinhos
df_simulacao$d <- NA
df_simulacao$kernel_type <- as.character(df_simulacao$kernel_type)
# funcao para paralelizar
func_llply <- function(i,data_frame=df_simulacao){
  df_temp <- data_frame
  sigma <- sigkernel(kernel = df_temp[i,"kernel_type"], 
                     p = df_temp[i,"k"], 
                     distance = df_temp[i,"dist_0"], 
                     density = df_temp[i,"DA"],
                     sigma.min=1e-6, 
                     sigma.max=1e6)$root  
}
registerDoMC(4)
replica.sim <- as.list(1:dim(df_simulacao)[1])
resultados <- llply(.data = replica.sim, .fun = func_llply, .parallel = TRUE)
df_simulacao$d <- unlist(resultados)
# padronização
df_simulacao %<>% mutate(k_perc = factor(k))
levels(df_simulacao$k_perc)[2] <- "0.10"
df_simulacao %>% names
df_simulacao %<>% select(SiteCode,p,k,DA,Ntotal,Stotal,txt.name,kernel_code,d,)
df_simulacao$txt.name <- gsub("dados_brutos","simulacao",df_simulacao$txt.name)

Figura 4 Distância média de dispersão em metros (d) ~ porcentagem de propagulos até a vizinhança imediata (k)

Simulação

A simulação coalescente foi construida em C++ e uma função em R foi utilizada para alimentar a simulação:

dinamica_coalescente <- function(U, S=0, N_simul, seed, disp_range, disp_kernel, landscape){
  # Runs coalescent simulations for a given heterogeneous landscape
  #
  # Parameters:
  # U: speciation rate
  # S: observed richness (integer) - used to fit the value of U, or set to
  #       0 (default) if that is not desired
  # N_simul: number of simulations
  # seed: seed of the RNG (an integer)
  # disp_range: width of the dispersal kernel
  # disp_kernel: an integer corresponding to the type of dispersal kernel. One of
  #               0: uniform
  #               1: normal
  #               2: Laplacian
  # landscape: either a filename containing the landscape data, or a
  #   bidimensional R array or matrix.
  #   TODO: describe the format of the input - trinary matrix)
  #
  # Returns:
  # r: an array of dimension N_simul x landscape dimensions, that is, each
  #   r[i.,] is a bidimensional array of the same shape as the landscape.
  #   Each site is labeled according to the identity of the species occupying
  #   that site.
  # U_est: estimated speciation rate. This is returned only if input parameter S > 0
  if (is.character(landscape)){
    l <- as.matrix(read.table(landscape))
    infile <- landscape
    land_dims <- dim(l)
  } else {
    land_dims <- dim(landscape)
    infile <- tempfile()
    # input file *must* be clean: no comments, headers or anything
    write.table(landscape, infile, col.names=F, row.names=F)
  }
  outfile <- tempfile()
  repeat {
    system(paste('./dinamica_coalescente', land_dims[1], land_dims[2], U, S, N_simul,
                 seed, disp_range, disp_kernel, infile, outfile))
    if (file.exists(outfile) || S == 0)
      break
    U <- U/10.
    print(paste("Decreasing value of U to", U))
    # set some lowest boundary here so simulations don't take forever
    if (U < 1e-20){
      print("Richness value too low, giving up...")
      return(NULL)
    }
  }
  r <- as.matrix(read.table(outfile))
  # transpose each grid, as output is written along lines but R reads it along columns
  # TODO: I thought I got it right, but it was wrong... please DO re-check
  #r <- aperm(r, c(1,3,2))
  
  # recover estimated speciation rate
  if (S > 0){
    out_con <- file(outfile)
    U_line <- strsplit(readLines(out_con, 2)[2], ' ')[[1]]
    close(out_con)
    U_est <- as.double(U_line[length(U_line)])
    return(list(r = r, U_est = U_est))
  }
  return(r)
}

Estimativa de U

Foram utilizadas 10 réplicas para estimar U:

# replicas
n_rep.U <- 10
func1 <- function(x,replicas=n_rep.U) {
  x$U <- NA
  x <- x[rep(1:dim(x)[1],each=replicas),]
}
df_referencia %<>% func1()
# simulacao
# valores de k
k_factor <- unique(df_referencia$k)
registerDoMC(n_cores)
for(a in 1:length(k_factor)){
  df_simU <- df_referencia %>% dplyr::filter(k == k_factor[a])
  op <- options(digits.secs=6)
  funcao_imigracao <- function(i,df_temp=df_simU){
    aviao <- list()
    aviao <- dinamica_coalescente(U = 1.25e-06, 
                                  S = df_temp[i,"Stotal"], 
                                  N_simul = 1, 
                                  seed = as.numeric(Sys.time()), 
                                  disp_range = df_temp[i,"d"], 
                                  disp_kernel = df_temp[i,"kernel_code"], 
                                  landscape = df_temp[i,"txt.name"])
    return(aviao$U_est)
  }
  replica.sim <- as.list(1:dim(df_simU)[1])
  sim.coal_U <- llply(.data = replica.sim, .fun = funcao_imigracao, .parallel = TRUE)
  df_simU[,"U"] <- unlist(sim.coal_U)
  write.csv(df_simU, 
            file=paste0("./U/","df_simU__k",k_factor[a],".csv"),row.names = FALSE)
}

Então calculou-se a média e variância por SiteCode e k (a maior variância foi na casa dos e-06). O U médio foi utilizado para parametrizar as predições de MNEE e MNEI (figura 6).

Figura 5 U_med estimado por k (% de prop até a vizinhança imediada) e d / L_plot (distância média de dispersão / largura da área amostral aproximada como um quadrado).

Estudo do padrão de U

Diferença entre pontos consecutivos

teste_acf <- df_resultados %>% filter(SiteCode==df_resultados$SiteCode[1] & MN == "EE") %>% .$U_med
hist()

SAD predita

MNEE

Para obter as SADs preditas por MNEE utilizei o seguinte código

f_simulacao <- function(i,df_=df_simulacao){
  X <- df_[i,]
  mat_sim <- dinamica_coalescente(U = X[,"U_med"], 
                                  S = 0, 
                                  N_simul = n_rep.SAD, 
                                  seed = as.numeric(Sys.time()), 
                                  disp_range = X[,"d"], 
                                  disp_kernel = X[,"kernel_code"], 
                                  landscape = X[,"txt.name"])
  l_SADs.preditas <- alply(mat_sim,1,function(Y) sort(as.integer(table(Y))) )
  file_name <- gsub("/home/danilo/Documentos/Doutorado/artigo_mestrado/1A.P_MOVER/simulacao/","",
                    X[,"txt.name"])
  file_name <- gsub(".txt","",file_name)
  path.file <- paste0(getwd(),"/SADs_preditas/",file_name,"__k",X[,"k"],".EE.", "rep_",1:length(l_SADs.preditas),".csv")
  for(j in 1:length(l_SADs.preditas)){
    write.csv(data.frame(SAD_predita = l_SADs.preditas[[j]]),
              file=path.file[j],
              row.names = FALSE)
  }
}
registerDoMC(4)
simulacao <- as.list(1:dim(df_simulacao)[1])
l_ply(simulacao,f_simulacao,.parallel = TRUE)

MNEI

Para MNEI:

# preparação dos dados
df_simulacao %<>% mutate(L_plot = 100*sqrt(Ntotal/DA),
                         m = d * ( 1 - exp(-L_plot*sqrt(2)/d) ) / (L_plot*sqrt(2)), 
                         m_ = m * p / (1 - (1-p) * m),
                         I = m_ * (Ntotal-1)/(1-m_),
                         J_M=p*DA*3600,
                         theta=(U_med*(J_M-1))/(1-U_med))
# simulação
f_simulacaoEI <- function(i,df_=df_simulacao){
  # i <- 1
  # n_rep.SAD <- 100
  df_name <- df_[i,]
  l_SADs <- replicate(n_rep.SAD,generate.ESF(theta = df_name$theta, I = df_name$I, J = df_name$Ntotal))
  file_name <- gsub("/home/danilo/Documentos/Doutorado/artigo_mestrado/1A.P_MOVER/simulacao/","",
                    df_name[,"txt.name"])
  file_name <- gsub(".txt","",file_name)
  path.file <- paste0(getwd(),"/SADs_preditas/",file_name,"__k",df_name[,"k"],".EI.", "rep_",1:length(l_SADs),".csv")
  for(j in 1:length(l_SADs)){
    write.csv(data.frame(SAD_predita = l_SADs[[j]]),
              file=path.file[j],
              row.names = FALSE)
  }
}
registerDoMC(4)
simulacao <- as.list(1:dim(df_simulacao)[1])
l_ply(simulacao,f_simulacaoEI,.parallel = TRUE)

Prévia Resultados

Calculo dos resultados

Para obter os resultados utilizei o seguinte código:

# agrupamento do df
df_auditoria %<>% group_by(SAD_obs.name) %>% nest
##### Resultados por Réplica #####
### rotina da função
df_auditoria$resultados <- vector("list",length = nrow(df_auditoria))
registerDoMC(4)
for(i in 1:nrow(df_auditoria)){
  # i <- 1
  df_ <- df_auditoria[i,]
  v_SAD.obs <- read.csv(df_$SAD_obs.name,header = TRUE,as.is = TRUE) %>% filter(species.correct != "Mortas") %>% 
    .$N %>% sort()
  df_predicao <- as.data.frame(df_$data[[1]])
  f_KSeS <- function(v_OBS = v_SAD.obs, path_SAD.MN){
    # path_SAD.MN <- df_predicao[1,"SAD_MN.name"]
    v_SAD.predita <- as.integer(read.csv(path_SAD.MN,header = TRUE,as.is = TRUE)$SAD_predita)
    a <- suppressWarnings(ks.test(x = v_OBS,
                                  y = v_SAD.predita))
    a <- data.frame(KS.D = a$statistic, KS.p = a$p.value)
    a$S_SAD.predita <- length(v_SAD.predita)
    a$S_SAD.obs <- length(v_SAD.obs) 
    return(a)
  }
  # f_KSeS(path_SAD.MN = df_predicao$SAD_MN.name[10])
  df_auditoria$resultados[[i]] <- adply(df_predicao,1, function(X) f_KSeS(path_SAD.MN = X$SAD_MN.name),.parallel = TRUE)
}
df_replicas <- df_auditoria %>% select(-data) %>% unnest(cols = c(resultados)) %>% as.data.frame()
##### Resultados #####
### rotina da função
alpha <- 0.05
df_resultados <- ddply(df_replicas,c("MN","k","SiteCode"),summarise,
                       D_mean = mean(KS.D), D_var = var(KS.D),
                       p.value_mean = mean(KS.p),p.value_var = var(KS.p),
                       S.MN_mean = mean(S_SAD.predita), S.MN_var = var(S_SAD.predita),
                       S.obs_mean=mean(S_SAD.obs),S.obs_var=var(S_SAD.obs),
                       n_SAD.N.ref = sum(KS.p>=alpha), n_SAD.ref = sum(KS.p<alpha),
                       .parallel = TRUE)

Auditoria 1: auto coerrência dos dados

Figura 6 Avaliação da autocoerrência dos métodos. geom_smooth(method=“loess”)

Padrões Gerais da congruência entre predito e observado

Figura 7 Padrões gerais para o número de SADs não refutadas (n_SAD.N.ref). geom_smooth(method=“loess”)

Figura 8 S_predita ~ S_obs. Linha vermelha: slope=1,intercept=0

Auditoria 2: Comparação com o conjunto de dados antigos

##                                   size isdir mode               mtime
## ./resultados/df_resultados.csv 1044347 FALSE  664 2019-11-14 07:55:10
##                                              ctime               atime
## ./resultados/df_resultados.csv 2019-11-14 07:55:10 2019-11-25 09:02:04
##                                 uid  gid  uname grname
## ./resultados/df_resultados.csv 1000 1000 danilo danilo

Figura 9 Diferença entre os parâmetros do novo conjunto de dados e do antigo.

Considerações
- o calculo de alguns parâmeros de MNEI foram modificado de 2017/2018 para o de 2019;
- os raster de paisagem e SAD observado foram atualizados;
- assim espera-se que: i) Ntotal_novo <= Ntotal_antigo & Stotal_novo <= Stotal_antigo & p_novo == | != p_antigo ii) d (distância média de dispersão para um determinado k), U e GOF_{MN:EE} sejam robustos;
Observado
parametros base
- Ntotal_novo >= Ntotal_antigo; Stotal_ e p_ idem;
parametros dos modelos
- d apresentou diff > 0, onde esperava que fosse uma distribuição simetrica centrada em zero (similar ao apresentado para U)
congruencia com o observado
- GOF_{EE}: a maior parte dos valores esta próximo de zero (50% dos dados) está próximo de 0, contudo há exemplos de valores onde houve total inversão nos valores (100 e -100)
- GOF_{EI}: a variância é superior a de EE, mas ainda a maior parte dos valores esta próximo de zero
- S.mean_{EE}:variância muito superior do que o esperado, uma vez que a regra é que EE ajuste corretamente S_observado
Perguntas
a) Quais os sítios que apresentam alteração em p, Ntotal e Stotal?
a.1) O erro encontrado nestes sítios pode ser avaliado para os 28 sítios adicionados?
b) Por que d apresentou diferença positiva?
c) Quais os sítios onde a congruência com o observado teve variação máxima? A hipótese de mudança de rotulos é plausível?

Auditoria dos dados

Sítios que apresentam alteração em p, Ntotal e Stotal

Figura 10 Paisagens onde houve diferença em p. 

• Sitios que podem corresponder a troca de rótulos por conta da diferença na proporção de habitat entre novo e antigo: MGuberl3, MGuberl5, MGuberl7, PEalia, PEmata2;
• Esses sítios apresentam alta proporção de habitat remanescente na paisagem e alta diff_p
• porém em alguns casos uma alteração na coordenada central pode levar à mudança de habitat observado: MGUberl1, BAlenc4, BAjuss e todos com diff_p -> 0
  
• pela minha avaliação visual a estimativa de p está correta, como havia avaliado anteriormente

SAD obs

Figura 11 Trabalhos que diferem em Ntotal e Stotal e proporção da abundância e riqueza das espécies indeterminadas (grep(" sp\.",species.correct)). Na primeira linha há os 75 pontos em comum (novo - antigo); na segunda linha há os 103 pontos do conjunto de dados novo. Novo: 2019; Antigo: 2017/2018

Sítios que apresentam alteração em GOF e S_mean

Figura 12 Comparação entre MN para as variáveis resposta (linhas 1 e 2) e para os parâmetros de dispersão (linha 3) por cenário de limitação de dispersão.

Congruência com o observado:GOF
- para alguns sítios independe o cenário de limitação à dispersão; porém para alguns há grande variação na resposta
- a maior parte dos sítios está próximo de zero para ambos os modelos, porém a variação de EI é superior.
Congruência com o observado:S.mean
- a riqueza média apresentou grande variação, inclusive em alguns sítios houve diminuição de 1000 spp do antigo para o novo - a diferença na riqueza média muda em função de k, o que não deveria ser observado em MNEE uma vez que independente de k ele estima corretamente a riqueza observada. Esse resultado é comum entre os dois conjuntos de dados (é intrinseco ao método). Assim essa mudança corrobora a hipótese de que pode existir uma troca de rótulos entre sítios. Parâmetros de Dispersão - a diferença em d pode ser explicada pelas mudanças aplicadas no conjunto de dados

Quais os sítios que devem ser avaliados

• Priorizar a avaliação da congruência para MNEE que deveria ser mais robusto
• Existe relação entre a diferença da congruência entre os conjuntos de dados e a proporção de espécies indeterminadas? Algum rótulo foi trocado?

Sítios para Avaliação Renato Lima